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Abstract 

Flutter instability in an infinite medium is a form of material instability corre- 
sponding to the occurrence of complex conjugate squares of the acceleration wave 
velocities. Although its occurrence is known to be possible in elastoplastic materials 
with nonassociative flow law and to correspond to some dynamically growing dis- 
turbance, its mechanical meaning has to date still eluded a precise interpretation. 
This is provided here by constructing the infinite-body, time-harmonic Green's 
function for the loading branch of an elastoplastic material in flutter conditions. 
Used as a perturbation, it reveals that flutter corresponds to a spatially blowing-up 
disturbance, exhibiting well-defined directional properties, determined by the wave 
directions for which the eigenvalues become complex conjugate. Flutter is shown to 
be connected to the formation of localized deformations, a dynamical phenomenon 
sharing geometrical similarities with the well-known mechanism of shear banding 
occurring under quasi-static loading. Flutter may occur much earlier than shear 
banding in a process of continued plastic deformation. 

KEYWORDS: Green's function; elastoplastic materials; nonassociative flow rule; mate- 
rial instability; granular materials. 
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1 Introduction 



Several micromechanisms acting at a microscale during deformation of granular and 
rock-like materials involve Coulomb friction. As a consequence, the flow rule becomes 
nonassociative and the phenomenological rate elastoplastic constitutive equations for 
these materials become unsymmetric. Due to this lack of symmetry, two squares of 
the propagation velocity of acceleration waves or, in other words, two eigenvalues of 
the acoustic tensor, may become a complex conjugate pair. That this situation might 
correspond to a form of material instability particularly relevant in granular material was 
clear since J.R. Rice (1977) coined for it the term 'flutter instability', but neither examples 
of constitutive equations displaying this instability nor a mechanical interpretation for it 
were given at that time. 

Consequently, research was initially focused on the determination of situations in 
which flutter was possible (see Bigoni, 2000; Loret et al., 2000 for reviews). In particu- 
lar, it was shown that flutter instability may occur more often than one might expect, 
not satisfying any hierarchical relation to other instabilities (such as for instance shear 
banding or second-order work negativity), possibly at an early stage of a hardening pro- 
cess and typically triggered by noncoaxiality (of the flow rule or induced by elastic or 
plastic anisotropy). However, the problem of finding a mechanical interpretation for the 
instability remained almost completely unexplored [with the exceptions of Bigoni and 
Willis (1994) and Sim5es (1997), the former considering a very simple problem setting 
and the latter providing some numerical tests]. This has been a major problem retarding 
further progress in research since, though generically believed to correspond to a dynam- 
ically growing disturbance, only the knowledge of the precise mechanical features of the 
instability can permit its identification for real materials. 

To shed light on this problem, a perturbative approach is developed in this article, 
following the methodology proposed by Bigoni and Capuani (2002; 2005) to investigate 
shear banding and other forms of material instabilities. In more detail, the analysis is 
limited in the present article to the loading branch^ of an elastoplastic constitutive oper- 
ator (taken from Bigoni and Petryk, 2002) embodying features typical of the behaviour 
of granular materials and capable of exhibiting flutter instability. An infinite body is 
considered made up of this material, homogeneously and quasi-statically deformed in 
two dimensions (plane strain or generalized plane stress). For this configuration a time- 
harmonic Green's function is found (in the way shown by Willis, 1991), which represents 
the first dynamic Green's function obtained for a nonsymmetric constitutive equation^. 

1 See Bigoni and Petryk (2002) for a discussion of this delicate assumption. 

2 A quasi-static Green's function for unsymmetric constitutive equation has been developed by 
Bertoldi et al. (2005), but this is unsuitable for flutter analyses, since this instability is essentially 
dynamic and thus remains unrevealed under the quasi-static assumption. In addition, Bertoldi et al. 
(2005) also derive boundary integral equations under the unsymmetric constitutive assumption, which 
are shown to possess certain typical features although not directly connected to the present discussion. 
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The Green's function is employed to form a pulsating dipole (two equal and opposite 
forces having a magnitude varying sinusoidally with time) to be used as a dynamic per- 
turbation revealing effects of flutter. 

Results demonstrate the following features of flutter instability that ma)| also occur in 
a material for which the tangent constitutive operator is positive definite (so that negative 
second-order work and shear bands are excluded at the considered stress level). 

• Differently from shear bands, becoming already evident when the boundary of the 
region of ellipticity is approached from its interior (Bigoni and Capuani, 2002; 
2005), flutter instability remains undetected while the eigenvalues of the acoustic 
tensor lie in the real range, appearing only after two real eigenvalues have coalesced 
and then become a complex conjugate pair; 

• flutter instability corresponds to a disturbance blowing-up in space from the per- 
turbing dipole and self- organizing along well-defined plane waves. 

• the normals to the above plane waves lie within the fan of directions corresponding 
to flutter and have been found to have an inclination remarkably different from 
that corresponding to shear bands, occurring later in the hardening process. 

It should be noted that the blow-up found in our solution will occur rapidly and 
nonlinearities neglected in our analysis (such as for instance the possibility of elastic un- 
loading and plastic reloading) may soon become important, possibly changing the overall 
mechanical response. Equally significant is the fact that the rate of growth increases 
with the frequency that is adopted. The governing equations of motion thus represent 
a problem that is dynamically ill-posed in the general transient case, unless the tangent 
moduli in fact display a frequency- dependence, such that the flutter effect reduces as 
frequency increases^]. However, our results suggest that flutter instability should induce 
a layering in an initially homogeneous material, inducing a localization of strain in a form 
somehow similar — though possibly occurring much earlier in a hardening process — to 
that pertaining to shear bands occurring in a dynamical context (Bigoni and Capuani, 
2005). Our hope is that this feature revealed by our results has now been made accessible 
to experimental investigation. 

1.1 Notation 

A standard, intrinsic notation is used throughout the paper (as for instance in Bigoni 
and Loret, 1999 and Bigoni, 2000), where vectors and second-order tensors are denoted 

3 More precisely, flutter instability has been shown by Bigoni and Loret (1999) to be unrelated to the 
occurrence of other instabilities such as loss of positive definiteness of second order work, loss of strong 
ellipticity and loss of ellipticity. 

4 Such a model was introduced by Bigoni and Willis (1994) in the context of a simple one-dimensional 
example. 
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by bold (the latter capital) letters. The scalar product between arbitrary tensors A and 
B is denoted by 

A - B — tr AB T = I • AB T , (1) 

where the usual symbols denoting the identity, the transpose, and the trace operator 
have been employed. In addition to the usual tensorial product between (vectors and) 
second-order tensors A and B 

(A ® B) [C\ — (B'C) A, (2) 
for every C, we will make use of the two tensorial products 

(Ap) [C] = A C + ° T B T , (A a B) [C] = ACB T , (3) 

so that I <g) J and J IE J become the symmetrizing and the identity fourth-order tensors, 
respectively. 



2 A simple constitutive model evidencing flutter in- 
stability 

We refer here to the model proposed by Bigoni and Petryk (2002) as a large strain 
version of that proposed by Bigoni and Loret (1999) [see also Bigoni (1995) and Bigoni 
and Zaccaria (1994)]. In particular, an objective symmetric flux, namely, the Oldroyd 
derivative of the Kirchhoff stress 

K= K — LK — KL T , (4) 

(where a dot over a symbol denotes material time derivative, L = FF^ 1 is the spatial 
velocity gradient and F the deformation gradient) is related to the Eulerian strain rate 

D = \{L + L T ), (5) 
through the piecewise-linear elastoplastic constitutive equation 

, E[D]-1(Q.E[£>])E[P] if/(K,/C) = 0, 
K= { H (6) 

E[D] if f(K,K.) < 0, 

where the symbol ( • ) denotes the Macaulay brackets operator (defined for every scalar a 
as (a) — (a+ \a\)/2), E is the elastic fourth-order tensor, / is the yield function in stress 
space depending on a collection K of internal variables (of arbitrary scalar or tensorial 
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nature); moreover, P and Q are the normals to the plastic potential and yield surface, 
respectively, and the plastic modulus H is related to the hardening modulus h through 

H = h + Q-E[P}. (7) 

In the present article, we will refer to the loading branch of eqn. ([6]), which is 

K = E[L) + LK+ K L - 4(E[P] ® E T [Q])[L], (8) 

where we have used the minor symmetries of E. Finally, introducing the first Piola- 
Kirchhoff stress 

S = KF T , (9) 



eqn. can be rewritten as 
where 



S = C[F], (10) 



(I ® F- r )E(I ® F~ T ) + I ® F- 1 S (11) 

- — (I El F- X )(E[P] g) E[Q])(J B F" T ). 
ii 



Note that the tangent constitutive operator C, eqn. ( fill) , possesses neither the minor nor 
the major symmetry, the latter except in the associative case, Q = P. 

2.1 Anisotropic elasticity 

Following Bigoni and Loret (1999) an anisotropic elastic law is assumed in the form 

E = \B®B + 2fj,B®_B, (12) 

where A and fi are material constants subject to the restrictions \x > 0, 3A + 2fi > 0, and 
B is a symmetric, positive definite second-order tensor, selected in the format 

B = b t b (g) b + b 2 (I -b®b), (13) 

where bi and b 2 are the eigenvalues of B, while the line spanned by the unit vector b and 
the plane perpendicular to it are the corresponding eigenspaces. Moreover, the material 
constants b\ and b 2 are assumed to depend on a single angular parameter b, restricted to 
the range ]0°,90°[ to meet the positive definiteness requirement of B, 

bi = v^cosfo, 6 2 = ^^si 11 ^ (14) 
so that the isotropic behaviour is recovered when b% — b 2 — 1, or b ~ 54.74°. 
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2.2 The acoustic tensor 

The acoustic tensor A ep (n) associated with the tangent constitutive operator C and the 
mass density p is defined by 

A ep (n)g = ^-C[g®n]n, (15) 

where n and g are the direction and amplitude of the propagating wave, respectively. 
Therefore, the acoustic tensor corresponding to C in eqn. (llip is 



A ep (n) = A e (n) - — (E[P]F~ T n ® E[Q]F~ T n) , (16) 
pH 

where A e (n) is the elastic acoustic tensor, defined as 

A e (n) = ^±^(BF- T n) <8> (BF~ T n) (17) 
P 

+ ^ \(F- T n) • (BF- T n)] B + -\n- (F^Sn)] I. 

p L J p L J 



Since C does not have the major symmetry, the acoustic tensor (|T6|) - (jT7|) is also not 
symmetric. 

2.3 Examples of nutter instability for plane problems 

The current configuration is assumed as reference, so that F = I and S = K = T, where 
T denotes the Cauchy stress. The plane problem is considered in which vector b and the 
propagation direction n lie in the plane spanned by ki and &2, two unit eigenvectors of 
K = T. Assuming the Drucker-Prager yield criterion, tensors P and Q take the form 

„ devT sinv , ^ , devT sinib , . . 

p = OTX ^Ti + vf 1 ' Q = ms,p \^T\ + (18) 

respectively, where devT = T — trT/3 and the angular parameters x an d ^ describe 
respectively the dilatancy and the pressure-sensitivity of the material. 

In the reference system {n, s,k 3 }, where s = k 3 x n, the acoustic tensor A ep (n) 
becomes 

/ A nn -^(n-q)(n-p) A e ns - -L{n ■ q)(s ■ p) \ 



A e ns -^(s.q)(n.p) At s --^(s.q)(s.p) 

pb 2 (n • Bri) + n • Tn 
P 



(19) 
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where 

q = E[Q]n = X(B ■ Q)Bn + 2pBQBn, 

(20) 

p = E[P]n = X(B • P)Bn + 2fiBPBn, 

and A e nnl A e ss , A ns are the in-plane components of the elastic acoustic tensor A e (n), 
namely 

A nn = X + 2fI (n • Bnf + -n • Tn, 

P P 

A e „ = ^-^(s • Bnf + -in • Bn)(s • Bs) + -n • Tn, (21) 
p p p 

A ns = ^^(n-Bn)(s-Bn). 

Note that the out-of-plane eigenvalue A33 in eqn. ( fl9l) corresponds to a wave with out- 
of-plane amplitude (g proportional to k^) and is assumed to remain strictly positive. 

From matrix ( TT9|) . we get the sum and the product of the two in-plane eigenvalues 
(squares of the acceleration waves propagation velocities) c\ and c 2 . corresponding to 
waves with in-plane amplitude (g lying in the plane spanned by kx and k 2 ), 

c\ + c\ = A e nn + A e ss —(fx - f 2 ), 

1 (22) 
c i c 2 = A e nn A e ss — (A e ns ) H 77(A e ns fz — A e ss fx + A nn f 2 ), 

pH 

where 

/1 = (n-q)(n-p), f 2 = -(a • q)(s • p), 

/a = (n-q)(s-p) + (s-qr)(n-p). 

A necessary and sufficient condition for the existence of complex conjugate eigenvalues 
af and is represented by the simultaneous fulfillment of the following three conditions 
(Bigoni and Loret, 1999) 

h = (A e nn - A e ss ) 2 [(fx + f 2 + 2e/ 3 ) 2 - (1 + 4e 2 )(/ 1 - f 2 f] > 0, 
/5 = (^„-^)(/ a + / 2 + 2e/ 3 )>0 > 



(24) 



h-^/u <p2R< h+vn 



where 



(^nn — ^L) 2 + 4(^ s ) 2 (Ann ~ A e ss ) 2 + A(A ns ) 2 

At. 



Ae Ae 

nn ss 



(25) 
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With reference to Fig. [TJ let Q a and 9 n be the angles of inclination of the direction of 
elastic anisotropy b and wave propagation normal n with respect to the stress principal 
axis fex- 




b 



Figure 1: Principal stress axes fci and £2, ax is of elastic symmetry b and propagation direction n, 
singled out by angles a and n , respectively. 

Dividing all quantities having the dimension of a stress in eqns. (|19p -(l24 p by fi, the 
parameters on which the condition of flutter depends are: 

• Elastic parameters: A///, strength of anisotropy b, and orientation of the axis of 
elastic symmetry with respect to the principal stress axis k±, namely 6 a . 

• Plastic parameters: plastic modulus H/fi, pressure sensitivity ip, and dilatancy x 
parameters. 

• Principal normalized deviatoric stress values: devTi/|devT|, devT 2 /|devT|, devT 3 /|devT 
However, these are not independent, so that given the form ( 1181) of P and Q, flutter 
depends on the angle 

/ devil devT 2 \ _ x ( [3 devT, \ 
L = sgn — — — + 2— — — cos \ —— — — 26 
& V|devT| |devT|y \ v V2|devT| y / v ; 

in the deviatoric plane, which is a 'modified Lode angle', defined for 61 G [— 7r, it] 
and in which sgn(O) = 1. 

It is possible to study flutter for all the propagation directions n while varying the 
plastic modulus H/fi and all remaining parameters in the above list are kept fixed, by 
use of inequalities (|24p . Therefore, the ranges in which flutter occurs can be plotted in 
the plane H/fi versus 9 n . Restricting the analysis to the infinitesimal theory, where the 
flux (jl]) is identified with T, analyses have been performed for simplicity with different 
values of the modified Lode parameter 6l = {60°, 30°, 0°, —30°, —60°}, as indicated in 

Fig. m 



s 




Figure 2: Stress directions in the deviatoric plane, denned by the modified Lode angle (j2l)]) . considered 
for flutter analysis. 



Results are reported in Figs. |3] and HI the latter giving more detail for four of the cases 
reported in the former figure. Different stress paths defined by the values of the modified 
Lode angle (I2"6"j) reported in Fig. [2] are considered for different anisotropy inclination Q a 
in Fig. [3] at given values of ip = 30° and x = 0°- m the graphs the closed contours 
denote regions where flutter occurs in the plane defined by the normalized critical plastic 
modulus H/fi and the inclination of propagation direction n . 

Four details of Fig. [3] are reported in Fig. HJ where A/// = 1, b — 80°, ip = 30°, and 
X = 0°, as in Fig. [3J The six regions in Fig. H] correspond to the four cases 9l = 0° and 
6 a = 15° (Case 1), B L = 6 a = 30° (Case 2), 6 L = and 9 a = 45° (Case 3), and 6 L = 
and (j = 60° (Case 4). 

With reference to the Cases 1,2,3 and 4, detailed in Fig. HI we note that the critical 
values of plastic modulus for loss of positive definiteness of the constitutive operator H^P 
and for loss of ellipticity Hf r permitting shear bands with normal inclined at 9 n E araj: 



Case 1: 




0.42, 




0.19, 


9nE = 


-28.0° 


Case 2: 


K D I» = 


1.22, 


H*/v = 


0.18, 


9nE = 


-16.4° 


Case 3: 


H cr D '/ ' ^ = 


1.03, 


H cr /n = 


0.74, 


OnE = 


-32.0° 


Case 4: 


H™/^ = 


1.84, 




1.57, 


QnE = 


-33.9° 



(27) 



so that in all cases flutter may initiate when the constitutive operator is positive definite 

5 Note that with 'ellipticity loss' we mean here the condition pertinent to the underlying quasi-static 
deformation. Moreover, due to anisotropy, only one shear band is found as first noticed by Bigoni et al. 
(2000). 
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(therefore at an early stage of a deformation process) and may extend in a region possibly 
involving loss of ellipticity. Note that thresholds ( 127)) have been graphically represented 
in Fig, HI where light grey regions correspond to regions where flutter may occur with 
the constitutive operator still positive definite, while in the dark grey regions ellipticity 
is lost (horizontal lines marking ellipticity loss are denoted with 'E (case i)', where i = 
1,..,4 stands for the number of the relevant Case). In the same figure, three black spots 
and a white spot (referred to Case 2) indicate the inclinations of shear bands at first loss 
of ellipticity. Note that the small flutter regions of Cases 3 and 4 are beyond the positive 
definiteness threshold, but still in the elliptic region. It may be important to remark that 

the initial inclinations of propagation normals for flutter and shear bands are 
unrelated and remarkably different. 

From the above analysis it can be deduced that the constitutive model allows one 
to approach flutter starting from a well-behaved state. Moreover, it may be interesting 
to note from Fig. H] that there are overlapping regions corresponding to different stress 
states (Cases 1 and 2). In these zones the flutter may have identical characteristics even 
if the stress state is different. 



2.4 Spectral analysis of the acoustic tensor 

The spectral analysis of the acoustic tensor is instrumental to the development of the 
Green's function that will be presented in the next Section. The analysis is restricted to 
the in-plane components of the acoustic tensor A ep 

A = A? x (k x ® k x ) + Af 2 {k x ® k 2 ) + A e 2 \{k 2 ® ki) + A c 2 p (k 2 <g> fc 2 ), (28) 

represented for later convenience in the principal stress basis k x , k 2 . The inverse of ( 12 8 j) 
can be written as 

A ^ = A e PA e P 1 A ep TTp ® *l) " ® **) (29) 

A XX A 22 /i 12 /i 21 

-A e 2 p (k 2 ® fej) + A e l{k 2 ® k 2 )} . 
Now, the eigenvalues of the acoustic tensor (128)) can be written in the form 

r 2 

2 | = /i n + ^ ±A ; A = y/(A? 1 -A%)* + 4A? 2 A1 



c 2 



1 21) 



(30) 



so that assuming non-defectiveness, the spectral representations of A and A 1 are 

A = c\{vi <g> Wi) + c\{v 2 ® w 2 ), (31) 
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Figure 3: Regions of flutter instability (occurring for internal points) in the H//i vs. 9 n plane, for 
the stress paths shown in Fig. [5] at various anisotropy inclinations 9 a . The following values of material 
parameters have been considered: X/fi = 1, b — 80°, ip — 30°, and x = 0°. 



and, assuming^] c\ ^ and c| ^ 0, 

= -gCUl <E> + ~2( v 2 ® *u 2 ), (32) 

C l C 2 

where {vi,^} and {tox, W2} are dual bases, thus satisfying Vi -Wj = Sij = 1,2), 

6 For A — > (coalescence of the eigenvalues), the tensor A becomes defective (except for the trivial 
case where A is isotropic) and each term in the spectral representation of A, and also of A -1 , blows up 
but A^ 1 continues to exist and to be defined correctly. Indeed a substitution of eqns. (|3"0T) and (j3"3")l or 
(|34j) into cqn. ((32]) leads to cqn. (l29l) . 
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Figure 4: Regions of nutter instability (occurring for internal points) in the H//J, vs. n plane, for 
A//i = 1, b = 80°, i\) = 30°, and X = 0°. Case 1: 6 L = 0° and 6 a = 15°. Case 2: 6 L = 30° and 
6 a = 30°. Case 3: as in case 1, but 6 a — 45°. Case 4: as in case 1, but 9 a = 60°. The regions of positive 
definitcness of the constitutive operator are marked in light grey, while (E) denotes loss of ellipticity into 
shear bands (regions shaded in dark grey) inclined at 6 n E(i), where i=l,..,4 denotes the relevant Case. 



composed of right, «j, and left, eigenvectors. This basis is given by 

«l — «1 H ?TTip K 2, «2 — K l H TTT^p K 2, 



2A A 2A A 



when 7^ 0, or by 

Z^i 2 l ^^21 

Wl = ~X kl + 2A fc2 ' ™ 2 = " A - * 1 + 2A * 2 ' 

when A% ^ 0. The case A% = Af 2 = is trivial. 



(33) 



(34) 
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3 The dynamic time-harmonic Green's function for 
general nonsymmetric constitutive equations 



An initial static homogeneous deformation of an infinite body is considered, satisfying 
equilibrium in terms of first Piola-Kirchhoff stress, namely, 

divS = 0, (35) 

and taken as the reference state in an updated Lagrangian formulation. A dynamic 
perturbation is superimposed upon this state, defined by an incremental displacement u 
satisfying the equations of incremental motion, written with reference to the constitutive 
equation ffTO]) in which dotted symbols are to be interpreted now as incremental quantities 
rather than rates. Thus 

CijkiUkjj + fi = pu ijtt , (36) 

where )t denotes material time derivative and fa and p are the incremental body forces 
and the mass density, respectively. 

Equations (!36|) look like ordinary elastodynamics, except that 

Cijki has neither the usual major ^ Ckuj nor the minor ^ ^ 
Cjiki symmetries. 

Note that tensor C^m can be identified (and will be in the examples) with that provided 
by eqn. ffTTj) . but can also be thought completely arbitrary in the following. To investigate 
the properties of eqn. fl3"6"|) . outside and inside the flutter region we follow the Bigoni and 
Capuani (2002; 2005) approach, based on the determination of the dynamic Green's 
function, sought for simplicity under the time-harmonic assumption 

Ui (x,t) = Ui (x)e- luJt , ffat) = fi(x)e-^\ (37) 

where u is the circular frequency and t and x denote time and space variables, respec- 
tively, so that the time dependence can be removed from eqn. (|36l) and consequently 

QjkiUkjj + p u 2 Ui + fi = 0. (38) 

The Green's tensor Gi P (x) is obtained by solving eqn. ( 138]) under the hypothesis fi = 
5i P 5(x), with 8(x) denoting the Dirac delta. We obtain 

Cij k iG kqt ij(x) + puj 2 G iq (x) + 8 iq 5(x) = 0. (39) 

In order to approach the flutter condition, we exploit the analysis of the acoustic 
tensor developed for the planar problem in Section 12. 3[ considering an infinite medium 
subject to plane strain (or generalized plane stress conditions), in which only four relevant 
components of the Green's function appear 

G iq = G iq (x 1 ,x 2 ), i,q = {1,2}, (40) 

and depend only on the two coordinates X\ and x^. 
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3.1 Radon transform 

The Green's function is determined employing a Radon transform technique [the alter- 
native approach employed by Bigoni and Capuani (2005) and based on a plane wave 
expansion is presented for completeness in Appendix A]. The Radon transform of a 
generic function f(x), x G R 2 is defined as 

K[f(x)] = f(p,n) = f{x)6(p-n-x)dx, p G R, n G R 2 (41) 

Jr. 2 



with the inverse 



m = ±f r/^-^s, (42) 
4vr 2 J\n\=iJ-oo (n-x-p) 

where a prime denotes partial differentiation in the following way 

/>,») = (43) 

In addition to the linearity, we will make use of the following properties of the Radon 
transform: 

• derivative transforms 

K [f,j(x)} = nj'ip, n), K [f^x)] = n in J"(p, n), (44) 

• transform of the two-dimensional Dirac delta function 

K[5(x)] = 5(p). (45) 

The Radon transform of eqn. is therefore 

CijkiriiUjGlgip, n) + pu 2 G iq (p, n) + 5 iq 5(p) = 0, (46) 

where 

Gl q {p,n) = —G kq {p,n). (47) 
Eqn. f|46|) can be rewritten in tensorial form as 

A(n)G"(p, n) + u 2 G(p, n) + = 0. (48) 

P 

Let us assume that A(n) has two non-null and distinct eigenvalues c 2 N and corre- 
sponding left and right eigenvectors w n, v n, (A = 1, 2), which can be used as dual basis 
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vectors, therefore satisfying v N -w M = 5nm, (N,M = 1,2). Employing the spectral 
representations of A(n) and I 



A(n) = 2J c 2 N v N <g) w N , I = v N <g) w N , 



N=l 



N=l 



in eqns. ( J48J) and representing the transformed Green's function as 



(49) 



G(p, n) = ^ <MP; n ) v N ® Wtf, 



AT=l 

where at is a (for the moment unknown) function of p and n, we get 

2 



AT=1 



P 



v N ® ty<v = 0, 



which is equivalent to the following uncoupled system of two equations, 



(50) 



(51) 



<>n + k 2 N <f>N + —$(p) = 0, N = l,2, 

P C N 



(52) 



where the wavenumber = uj/cn has been introduced. Since we have chosen the 
harmonic time dependence to be of the form e~ luJt , the outgoing wave solution of (152]) in 
the p coordinate is: 



(f) N (p, n) 



ik N \p\ 



2 • 



2pik N c N 



so that 



and 



- e ik N \p\ 



G(p,n) = -Y,^—TV N ®w N . 



N-- 
2 



#(p,n) = -£ 



J x 2pik N c N 
sgn(p)e ifcjv ' p ' 



N=l 



2pc 



2 -V N ($W N . 

N 



(53) 



(54) 



(55) 



The antitransform of equation ( 1541) leads to 



G(x) 



An 2 



^2 



N=l 



n 



+ °° sgnQoW^I , , 
Ytz — : ~ v n ® w Nap as. 



ii-00 2pc N (n-x-p) 



(56) 



The integral in the variable p can be evaluated in the way shown in Appendix B, so that, 
employing the cosine and sine integral functions 



Ci(z) 



z cost 



dt, |argz|<7r and Si(z) 



sint 



dt. 



(57) 
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the Green's function can be finally written in the form 

G(x) = --^ V I [2cos(k N n ■ x)Ci(k N \n ■ x\) (58) 

+2sm(kj^n • x)Si(k]yn • x) — iTccos^kj^n • x)\ — ^ 2 — —ds. 

P C N 

We introduce polar coordinates so that the position vector x has modulus r = |ai| and is 
inclined at angle 9 to the Xi-axis. Taking the unit vector n inclined at a + 9 with respect 
to the xi-axis (so that a is the angle between x and n) and noting that cos(-) Ci(-) and 
sin(-) Si(-) are even functions, we can re-write eqn. ( 1581) as 

G(x) = — — — ■ } / [2cos(r/cAr|coso;|)Ci(r/c7v|cosa|) (59) 

N=l J0 

+2sin(r/cAr|cosa;|)Si(r/cAr|cosa;|) — 27rcos(rA;Ar|cos a\)] — ® 2 — —da, 

where k^, v N , w N and c 2 N depend on a + 9. 

The acoustic tensor is a periodic function of a with period 7r since 

pA ik (n) = CuklTll + (Cufc2 + Q2fcl)^1^2 + Q2fc2«2, (60) 

where n\ = cos(a + 9) and n 2 = sin(a + 9), and also cat, k^, «jv, and iuat are periodic 
functions of a with the same period. It follows that the integrand in eqn. fl59|) is 7r- 
periodic. Therefore, 

the two-dimensional, time- harmonic Green's function corresponding to a generic, 
completely non- symmetric constitutive fourth-order tensor, relating the incre- 
ment of the first Piola-Kirchhoff stress to the deformation gradient increment, 
eqn. fW\) , can be written in the form 



1 y 2 r 

G(x) = -y / [2cos(r/cAr|cosa;|)Ci(r/cAr|cosa|) (61) 

An 1 In 

N=l Jl> 

+2sin(rfcAr|cosa|)Si(rA;jv|cosa|) — z7rcos(r/cAr|cosa:|)] — < ^ > 2 — —da, 

P C N 

where k^ = u/cn and c% are the eigenvalues of the acoustic tensor A, eqn. (|3T1) with 
corresponding left and right eigenvectors and v^, all quantities depending on n, 
which means on a + 9. 

It can be noted that the integrand in eqn. flSTj) displays a logarithmic singularity at 
r = and a = ir/2, since (Lebedev, 1965) 

f z 1 — cos t 

Ci(z) = 7 + log z — / dt, |argz|<7r, (62) 

Jo t 

where 7 is Euler's constant. 
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4 A dynamical interpretation of flutter instability 



The dynamical interpretation of nutter instability will be achieved following the approach 
introduced by Bigoni and Capuani (2002; 2005), so that the Green's function is employed 
to provide a dynamical perturbation to be superimposed upon a given state of equilibrium 
of a homogeneously deformed material. Several plots of Green's tensor components will 
be presented, so that a preliminary normalization of the Green's tensor and a study of the 
involved non-dimensional parameters becomes instrumental. In particular, introducing 
an arbitrary characteristic length a and consequently the dimensionless spatial variable 
x = x/a, making use of the property 

5 (ax) = -z6(x), x e R 2 , (63) 
a 2 

eqn. ( 1391) can be rewritten as 

Cijk fP^ + u 2 G tq (x) + S ig S(x) = 0, xeR 2 (64) 



where 



n v v 

Thus, a dimensionless version of the Green's tensor fl58l) reads 



Qjki = — UJ = a \rr. UJ - ( 65 ) 



1 2 f 

G(x) = 2/ / [2cos(k N n • x)Ci(k N \n • x\) (66) 

+2sm.{k^n • x)Si{kj^n • x) — iircos(kN n ' x)] — % — ~ 



where 



k N = ak N = —, c N = J-c N , (67) 
c N V A* 



so that c N are the eigenvalues of the dimensionless acoustic tensor A = pA/fj,. 

4.1 Effects of flutter instability on Green's tensor 

The behaviour of the Green's function, eqn. flBT]) . is briefly analyzed here, outside and 
inside the flutter region. As a reference, we consider Case 3 shown in Fig. HI in which the 
material is subject to the radial stress path corresponding to 6l = in Fig. |2] and the 
direction of the axis of elastic symmetry is taken inclined at 8 a = 45° with respect to the 
principal stress direction fei. The employed material parameters are A/// = 1, b — 80°, 
ijj = 30°, and x — 0°. The dimensionless Green's tensor components have been computed 
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for u) — 1 and for several values of the plastic modulus H/fi, including the values 3.53, 
and 1.5. These correspond, respectively, to situations near and inside the flutter region 
(see Fig. H|), but still in a situation where the constitutive operator is positive definite. 
The values of the components are plotted in Fig. [5] as functions of the distance from the 
singularity along a radial line inclined at —45° with respect to the x\ axis, normalized 
through division by a. The real (imaginary) parts of the Green's function components are 
plotted left (right) in the figure, the plots having been obtained starting from x\ = 1/10 
to exclude the singularity (in the real components of the Green's tensor). 

Commenting on the results, first, we note from the figure that the Green's tensor is 
not symmetric (since the acoustic tensor is not), so that Gu ^ G^i- 

Second, results referring to values of plastic modulus H//J, higher than 3.53 and up 
to 7, not reported here for conciseness, produce curves practically coincident to those 
pertaining to H/fi = 3.53; we can therefore conclude that there is not much difference 
between the situations in which the material is far from and very near to the flutter 
region. This feature has been confirmed by us with several calculations (not reported 
here) and distinguishes flutter from shear banding, the latter becoming already visible 
when the condition of loss of ellipticity is approached from the interior of the elliptic 
range (Bigoni and Capuani, 2002; 2005). 

Third, a blow-up of the solution with the space variable, clearly visible in all com- 
ponents of the Green's tensor is the characteristic feature of instability inside the flutter 
region, H/fi = 1.5. This blow-up is similar to that evidenced by Bigoni and Willis (1994), 
but in a constitutive setting including viscosity, which is now absent. 

It becomes evident that further exploration of flutter instability requires plotting of 
incremental displacement maps. These are obtained below employing a perturbation in 
the form of a pulsating dipole. 

4.2 Effects of flutter instability revealed by a perturbing dipole 

The singular solution previously obtained, eqn. ( loTj) . can be used to analyze the effects of 
a perturbation superimposed upon a given homogeneous deformation of an infinite body. 
We follow here Bigoni and Capuani (2005) considering the simplest self-equilibrated per- 
turbation in terms of a dipole: two equal and opposite pulsating forces of unit amplitude, 
taken at a distance 2a apart, along a line inclined at /3 = 45° with respect to the Xi-axis, 
see Fig. El 

For this loading system, the level sets of the real part (left in the figures) and the 
imaginary part (right in the figures) of the components U\ (first and third parts from 
the top of the figure) and u% (second and fourth parts from the top of the figure) of 
incremental displacements have been computed and plotted in Figs. I7HT21 The two 
upper parts of all the figures refer to a situation far from flutter instability, whereas the 
two lower parts refer to a situation of flutter, well inside the region of instability. 
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Figure 5: Dimensionless Green's tensor components (real part left, imaginary part right in the figure) 
along a radial line inclined at —45° with respect to xi-axis, for Case 3 of Fig. @] and oj = 1. Two values 
of the plastic modulus H / '/i = {3.53, 1.5} are considered, corresponding, respectively, to situations near 
and inside the flutter region. The blow-up of all components of the Green's tensor is evident in the 
flutter region, H/ji = 1.5. 
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Figure 6: Geometry of the time-harmonic pulsating perturbing dipole. 



The following parameters have been selected to be equal for all figures: 

A//i=l, 6 = 80°, ^ = 30°, x = 0°. 

Moreover, Figs. [7HT01 refer to the same nondimensional frequency parameter u> — 1, 
whereas the effect of frequency is explored in Figs. [11] and [121 pertaining respectively to 
u = 2 and 1/2 and corresponding to the same parameters employed in Fig. [81 All compo- 
nents of incremental displacements have been plotted for the nondimensional coordinates 
Xi/a and x%ja ranging between —25 and 25, with the exception of Fig. [101 where this 
range has been extended to —50 and 50 to help visualization of the blowing-up typical 
of flutter. 

The differences between Figs. ITrfTOl lie in the choice of different stress states expressed 
in terms of 9i and anisotropy direction 9 a . In particular: 

• Fig. [prefers to H/ \i = 3 (two upper parts), H/fi = 0.32 (two lower parts) and to 
Case 1 of Fig. gj where 9 L = 0° and 9 a = 15°; 

• Fig. prefers to H/fi = 2 (two upper parts), H/fx = 0.25 (two lower parts) and to 
Case 2 of Fig. HI where 6 L = 30° and 9 a = 30°; 

• Fig. prefers to H/ji = 4 (two upper parts), H/jj, = 1.5 (two lower parts) and to 
Case 3 of Fig. HI where 6 L = 0° and 9 a = 45°; 

• Fig. [TU1 refers to H/ ji = 4 (two upper parts), H//i = 1.9 (two lower parts) and to 
Case 4 of Fig. H where 9 L = 0° and 9 a = 60°. 

Note that the values of the plastic modulus selected for the examples are all higher 
than the critical values for loss of ellipticit;y0 [see the values listed in (I27p ). so that shear 

7 More precisely, all the considered plastic moduli are higher than the critical values for loss of strong 
ellipticity (Bigoni, 2000). 
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Figure 7: Level sets of the real (left) and imaginary (right) parts of the components of incremental 
displacements {u\ first and third parts from the top, u-i second and fourth parts) for a dipole inclined 
at /3 — 45°, far from (upper two parts, H/n = 3) and inside (lower two parts, H/fi = 0.32) the flutter 
region. Results pertain to Case 1 of Fig. |H for ui — 1. Note the system of blowing-up, parallel waves 
revealing the effect of flutter. 
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Figure 8: Level sets of the real (left) and imaginary (right) parts of (the modulus of) incremental 
displacements for a dipole inclined at /? = 45°, far from (upper part, H//.1 = 2) and inside (lower part, 
H/fi, = 0.25) the flutter region. Results pertain to Case 2 of Fig.0] for uj — 1 . 
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bands are excluded. However, all the values of H/fi corresponding to situations far 
from nutter and the two values 1.5 and 1.9 lie in the zone of positive definiteness of the 
constitutive operator, while the two values 0.25 and 0.32 have been selected outside this 
region [see the values listed in (127)) ]. 

It can be observed from the upper parts of Figs. ITHTOl (referring to a non-flutter 
situation) that the displacement maps are typical of an anisotropic material, since 45°- 
symmetry is not in evidence. Moreover, decay of the solution is appreciable, when the 
distance from the dipole increases. Now, considering the lower parts of the figures, the 
effects of flutter instability become self-evident. In particular, we may observe a growth 
of the solution in space, which tends to degenerate into a system of blowing-up, parallel 
plane waves. Results not reported here for brevity demonstrate that: 

the inclination of the blowing-up plane waves is almost independent of the 
dipole inclination ( angle (3 in Fig. so that it has to be considered a charac- 
teristic of the material, related to the particular stress state and constitutive 
features. We have observed that the inclination of the blowing-up waves cor- 
responds to a value in the middle of the inclination fan of flutter (see Fig. \3^i. 

In particular, the inclinations of the plane waves at a sufficient distance from the 
dipole, are different in Figs. [TUTUl but correspond to the mean value of flutter direction fan 
visible in Fig. H] at the analyzed Hj \x values. On the other hand, the same inclinations are 
found for figures Figs. [S] and [TT] and H21 since these cases differ only in the nondimensional 
frequency parameter ui, which influences only the spacing of the blowing-up waves. 

As far as the effects of varying the nondimensional frequency parameter u> are con- 
cerned (see Figs. Hi] and H2J referring to the same material parameters as in Fig. [HI but 
with Co = {1, 2, 1/2}), we see that an increase in the frequency yields a narrowing of the 
distance between blowing-up plane waves. Moreover, increase in frequency gives rise to 
the 'shadowing' effect already noted by Bigoni and Capuani (2005) for shear bands. 

Compared to the shear bands analyzed by Bigoni and Capuani (2002; 2005), we may 
observe that these are already revealed when the boundary of the region of ellipticity 
is approached from the inside, while flutter remains undetected. Beside this difference, 
there are however many similarities between the two phenomena: first of all, shear bands 
tend to blow-up in space as the boundary of instability is approached, and extend from 
a perturbation to infinity, outside the elliptic range. Second, shear bands also tend 
to degenerate into families of plane waves parallel to a specific direction. Third, the 
signals tend to focus along well defined patterns, both for shear bands and for flutter. 
Note however, that flutter instability may occur much earlier than shear banding in a 
deformation process; moreover, waves near the loss of ellipticity threshold tend to blow- 
up along the shear bands but, in contrast to flutter, they tend to decay in the parallel 
direction. 

As a conclusion, we remark that flutter instability yields a self-organization of dy- 
namic disturbances along well-defined and blowing-up parallel waves, having inclinations 
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corresponding to the mean value of the inclinations for which flutter is possible at the 
considered constitutive setting and stress state. 

From the mechanical point of view, our results suggest that flutter yields a 'layering' 
of deformation patterns, with an inclination corresponding to the flutter direction, a 
spacing related to the frequency of the perturbing agency, and possibly occurring early 
in a plastic deformation process. 

5 Conclusions 

Following the approach to material instabilities proposed by Bigoni and Capuani (2002; 
2005), flutter instability in a continuous elastoplastic medium has been investigated, by 
finding the dynamic, time-harmonic Green's function for the loading branch of a fully 
unsymmetric tangent constitutive operator, embodying features typical of the behaviour 
of granular materials. For this material, flutter instability may occur when the constitu- 
tive operator is positive definite (so that the solution of the rate infinitesimal problem is 
unique and shear bands are excluded), while two eigenvalues of the acoustic tensor are 
complex conjugate. Our results provide the first interpretation of flutter instability, which 
is shown to correspond to a dynamical instability growing in space and self-organizing 
into plane waves with normals lying in the fan corresponding to the complex eigenvalues 
of the acoustic tensor and yielding a sort of 'layering' of unstable deformation patterns, 
showing some similarity to shear band instability. The rate of growth of the solutions 
displayed here increases with the frequency that is assumed. This demonstrates dynam- 
ical ill-posedness of the governing equations of motion in the general transient case and 
implies a need that is physical as well as mathematical for the admission of some appro- 
priate rate-dependence into the constitutive model, to remove the flutter effect at high 
frequencies. Although no such mechanism is built into the present analysis (the tangent 
moduli would become functions of u but this is in any case fixed), and other mechanisms 
not accounted for (such as for instance the possibility of elastic unloading and material 
viscosity) may change some of our conclusions, we believe that the emergence of the 
layered structures that we have found may find future experimental validation. 
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Figure 9: Level sets of the real (left) and imaginary (right) parts of (the modulus of) incremental 
displacements for a dipole inclined at ft = 45°, far from (upper part, H/ ' \i = 4) and inside (lower part, 
H/fj, = 1.5) the flutter region. Results pertain to Case 3 of Fig. SI for Q = 1 . 
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Figure 10: Level sets of the real (left) and imaginary (right) parts of (the modulus of) incremental 
displacements for a dipole inclined at j3 = 45°, far from (upper part, H/ fx = 4) and inside (lower part, 
H/ /i = 1.9) the flutter region. Results pertain to Case 4 of Fig. HJ for 6D = 1 . 
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Figure 12: As for Pig. EJ but with Q = 1/2 . 
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APPENDIX A. Green's function obtained via plane wave expansion. 

The Green's function (I58I) is obtained here for completeness using the plane wave 
expansion technique employed by Bigoni and Capuani (2005). The plane wave expansion 
of the 5 function and of the Green's tensor G(x) are, respectively, 



5(x) = -— -da, G(x) = -— G(n-x)ds } (A.l) 

4vr 2 J m=1 (n • xf 4vr 2 J m=l 

where n is a unit vector, so that the plane wave expansion of eqn. (|39|) leads to 

C ijkl n jni Gl q (Z) +pu 2 G iq (Z) + |i = 0, (A.2) 

where £ = n ■ a:. In this equation the acoustic tensor can be easily recognized, = 
CijkiTijni, so that we get 

A(n)G"(0 + co 2 G(0 + ^1 = 0. (A.3) 
Writing now the analogue of the representation ( l50l) . namely, 

2 



G(0 = <M0«Jv ® iwjv, (A.4) 

AT=1 

we transform eqn. ( 1A.3I) into the analogue of eqn. ( 15"TT) 

y] (c 2 N (j)" N + W 2 07V + -TJ ) «JV ® «>JV = 0, 
2V=1 v Ps / 



(A.5) 



which is equivalent to the following uncoupled system of two equations, analogous to 
eqns. ([52]), 

<Pn + k 2 N (f> N + \\- 2 =0, N = 1, 2, (A.6) 

where k^ = oj/cn. 

The sole physically meaningful solution of the ordinary differential equation (IA.6I) is 
obtained by imposing the radiation condition, stating that the solution should include 
only outgoing waves. Since the harmonic time dependence has been selected in the form 
e~ tujt , the outgoing wave solution of (I A . 6 [) in the £ coordinate is: 

<M0 = 77^T [2Ci(k N \^\)cos(k N (A.7) 
z,pc N 

+2Si(/cAr£)sin(fcjv£) — 27tcos(/cat£)] . 
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Finally, a chain of substitutions, of eqn. flA.7[) into eqn. ()A.4j) and finally into eqn. 
(jA.ip ?. leads to the Green's function in the form (|58|) . 



APPENDIX B. Evaluation of the integral in the variable p in eqn. ( 156ft . 

The integral in the variable p appearing in eqn. ( )56|) can be evaluated splitting the 
domain as follows 

sgn W e d p=- e _ dp+ | dpj (R1) 



z-p J-oo i-p Jo t-v 

so that we can treat the two integrals separately, namely 

p -ik N p r+oo iq 

dp = -e~ ik "Z / — dq, (B.2) 

oo S — P Jk N £_ Q 

where we have made the substitution q = &jv(£ — p), and 

+00 ik N p r+oo iq 

dp = -e ikN ^ / — dq, (B.3) 

i-p J-k N z q 

where we have made the substitution q = fcftr(p — £). The two expressions (IB. 21) and (1B.3[) 
are used to get eqn. (I58"|) . 
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